A pseudo-spectral approach to inverse problems in interface dynamics 
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An improved scheme for computing coupling parameters of 
the Kardar-Parisi-Zhang equation from a collection of succes- 
sive interface profiles, is presented. The approach hinges on a 
spectral representation of this equation. An appropriate dis- 
cretization based on a Fourier representation, is discussed as a 
by-product of the above scheme. Our method is first tested on 
profiles generated by a one-dimensional Kardar-Parisi-Zhang 
equation where it is shown to reproduce the input parame- 
ters very accurately. When applied to microscopic models of 
growth, it provides the values of the coupling parameters as- 
sociated with the corresponding continuum equations. This 
technique favorably compares with previous methods based 
on real space schemes. 

64.60.Ht,05.40.+j,05.70.Ln 



I. INTRODUCTION 

In most complex systems, it is often difficult to di- 
rectly relate microscopic interactions to the dynamics of 
coarse-grained (mesoscale or large scale) spatial struc- 
tures. In this context, nonlinear inverse methods |l]] 
which infer equations governing a system from experi- 
mental observations of its successive time evolution, may 
prove to be more efficient then direct methods. When 
only experimental data coupled with the hypothesis of 
an underlying determinism are used, the identification is 
purely non-parametric. Such methods are extensively ex- 
ploited for instance in biological and economical systems 
where predictions are typically relying not on basic mech- 
anisms, but directly extrapolated from time series using 
various procedures (neural networks, nearest-neighbors 
algorithms, etc) M . On the other hand if a parametrized 
phenomenological equation, usually derived from a com- 
bination of general symmetry considerations and heuris- 
tic physical arguments, is assumed from the outset, the 
inverse approach consists in finding an optimal set of pa- 
rameters. In this work, we focus on this latter situation 
in the framework of interface growth dynamics. 

The phenomenon of interface growth covers many tech- 
nological applications ranging from epitaxial deposition 
to bacterial growth and fluid motion in porous me- 
dia |^-^|. It is known that the large scale dynamics of 
such systems may be conveniently address in terms of 
continuum stochastic differential equations. The Kardar- 
Parisi-Zhang (KPZ) equation constitutes a paradig- 
matic universality class which is believed to aggregate 



a large portion of observed interface dynamics. This 
Langevin type equation possesses a non-linear term 
which accounts for the interface local normal growth ab- 
sent in its linear counterpart the Edward- Wilkinson equa- 
tion (EW) H . Despite a major effort both from the nu- 
merical and analytical point of view, a complete char- 
acterization of the KPZ properties is nevertheless still 
lacking. From the numerical viewpoint, finite-difference 
schemes have been widely exploited to approximate the 
continuum KPZ dynamics. Unfortunately, naive dis- 
cretizations of the spatial derivatives may result into be- 
haviors inconsistent with known properties of the con- 
tinuum equation. For instance, the usual symmetric 
second-order finite difference scheme for the non-linear 
term violates an important symmetry of the continuum 
one-dimensional theory |lC[ |. An appropriate modifica- 
tion in the framework of finite difference approximations 
may heal this drawback 0,0, but it is nonetheless un- 
able to properly display coar se-g rained properties of the 
original continuum equation [pL 2|| . Specifically it does not 
preserve the correct functional form of the coarse-grained 
equilibrium distribution, a basic feature that one expects 
from the renormalization group point of view (a similar 
feature was already observed in Ref |l3| ). 

In the present work, we follow a different route by 
introducing a numerical approach which preserves more 
features of the original continuum KPZ equation. From 
the outset, our method is based on a spectral rather than 
a finite difference scheme. Spectral methods are widely 
used in fluid mechanics |bf | and their accuracy and relia- 
bility compared to those of finite difference schemes have 
been tested, over the years, on deterministic equations. 
Either in the context of fluctuating interfaces, it has al- 
ready been used |l^,[l6| to analyze the evolution of a de- 
terministic equation, namely the Kuramoto-Shivashinsky 
equation. However, it has not been directly implemented 
on the stochastic partial differential equation |17|. In 
this paper, we first extend spectral methods to stochastic 
equations, using the KPZ equation as a testbench, then 
devise a reconstruction procedure for the inverse problem 
for the KPZ, and finally show that it is unaffected by the 
deficiencies of real space approximations. 

The first reconstruction of the stochastic KPZ dynam- 
ics from experimental surface has been performed by 
Lam and Sander pq . These authors used a standard 
finite-difference scheme to approximate the dynamics and 
based their reconstruction on a least-squares algorithm. 
Experimental dynamics was obtained by simulating var- 
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ious microscopic models. This strategy, however, seems 
to be limited for two reasons. First it is based on a 
problematic numerical approximations of the KPZ equa- 
tion. Second the least-squares algorithm, originally tai- 
lored for deterministic equations was directly transpose 
to stochastic or Langevin type equations |l8[ . In fact, 
while the performance of this method has been tested in 
the presence of measurement noise - noise which affects 
experimental observations but does not change determin- 
istic trajectories - in the presence of dynamical noise (e.g. 
Langevin dynamics) this scheme can be far more unsuc- 
cessful. Moreover, the least-squares algorithm of Ref. 
is very much dependent on the assumption of small sam- 
pling time of observations a fact that should be checked 
a posteriori (see Ref. [12] for more details). 

In a previous paper (13), an alternative reconstruc- 
tion algorithm which does not suffer of the aforemen- 
tioned problems, was introduced. The basic strategy 
amounts to apply the least-squares procedure at the level 
of correlation functions rather than directly to the inter- 
face stochastic variables. This method, albeit successful, 
could not correctly account for the surface coarse grained 
properties because of an intrinsic deficiency of spatial dis- 
cretizations. Motivated by this feature, we shall remove 
this drawback by introducing a spectral representation of 
the Langevin equation which is not plagued by discretiza- 
tion problems. Then we shall proceed to reformulate and 
test the approach of Ref [12] for the inverse KPZ prob- 
lem, in spectral space. 

The paper is roughly divided in two parts and orga- 
nized as follows. The first part concerns the KPZ equa- 
tion itself and its approximations. In Section ||, the con- 
tinuum KPZ equation in (1 + l)-dimcnsion is rapidly re- 
called along with some previous spatial discretizations 
used to numerically mimic it. In the same Section, a 
novel discretization in Fourier space is then proposed and 
shown to improve over the real space approach from a 
purely theoretical viewpoint. The deterministic evolution 
equation for correlations functions in spectral space, are 
then derived in Section III. Section [V is devoted to the 



numerical procedure and, in particular, to the treatment, 
via a pseudo-spectral method, of the nonlinear term of 
the KPZ equation. Numerical results are then provided 
to compare the performance of the various discretizations 
to the corresponding analytical results obtained for the 
continuum KPZ equation. The second part of the pa- 
per is devoted to the inverse method for KPZ dynamics, 
where we discuss how one can reconstruct the dynamics 
from the knowledge of interface profiles. This technique, 
based on equations derived in Section [H]|, is described in 
Section |v|. It is then tested in Section VI against data 
produced by synthetic interface profiles generated by a 
(1 + l)-dimensional KPZ equation with known coupling 
parameters. The test is performed even in the presence 
of coarse-graining and provides the renormalization prop- 
erties of the KPZ coupling parameters. Finally, in Sec- 
tion VII , the reconstruction method is applied to various 



II. SPECTRAL DISCRETIZATION OF THE KPZ 
EQUATION 

We consider a one-dimensional surface profile of width 
L. The surface can be either experimentally or numer- 
ically generated, and we assume it to be periodic with 
period L. At a mesoscopic scale, it can be described 
within a continuum (hydrodynamic) approximation by 
a variable h(x,t) (with < x < L) which satisfies a 
Langevin equation A paradigmatic example is given 
by the Kardar-Parisi-Zhang (KPZ) equation Q: 



d t h = c + vdlh+^{d x hf + r,, 



(1) 



where ri(x, t) is a noise with zero average and (^-correlated 
both in time and space as: 



{rj(x, t)r)(x', t')) = 2DS(x - x') 5(t - t'). 



(2) 



In Eq. (|), symbol (...) means that an ensemble aver- 
age over the noise is performed and c, v, A, and D are 
coupling parameters. In writing Eqs.(|l) and (||), an ap- 
propriate regularization is always tacit y assumed. This 
amounts to define a minimal length scale a for h(x,t) 
and to assume that the surface variable h(x,t) satisfies 
a KPZ regularized equation with a given spatial cut-off 
depending on scale a. For instance, this can be done by 
considering a discretization version of Eqns. (0) and (|) 
with a as a lattice constant : 



dh 
It 



_ =c+ ^ Fm+ A. m h}+\ - 



or alternatively 



dt 



2a 2 



(3) 



(4) 



where hj stands for the value h(xj,t) of the periodic 
smoothed surface at Xj = ja, j = 1,...,N where 
N = L/a and the quantities 9i(t) are uncorrelated white 
noise functions 



{e i {t)9 J (t')) = 25 ij 8{t-t') 



(5) 



microscopic models at various coarse-grained scales. 



In Eqs. (g) and (Q), one usually approximate the dis- 
cretized Laplacian with the second-order finite difference 
term 

F?[h] = h i+1 + h i -i-2h i , (6) 
and the non-linear term with 

^ A N = ^+i-^-i] 2 , (7) 

or 

K X [h} =l[(h i+ i - hi) 2 + (hi - h l+1 f (8) 
+ (h i+ i - hi)(hi - hi-x)] . 
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As previously discussed |1Q-|12|, only choice given by 
Eq. ([|) guarantees that at least some properties of the 
correct equilibrium distribution are retrieved. Both rep- 
resentations Eqs. (Q) and (^) are however problematic 
at the coarse-grained level |12] . 

Here we show how one can achieve an alternative and 
always correct discretization of the continuum KPZ equa- 
tion by a procedure that directly applies in momentum 
space. 

We first expand the periodic continuum field h(x,t) in 
Fourier modes as 
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where the Fourier component 

oL/2 



dx h(x, t) e- iq " 



(9) 



(10) 



-L/2 



is associated with wavenumber q n = 2im/L. Since h(x, t) 
is real, this imposes h* n — h- qn or alternatively, if h Qn = 



<-q„ 



-i/3 qn is separated in real and imaginary part, S_ 9n = 



-j3 q . A similar expansion for the noise term 



r](x, t) leads to Fourier components rj qn with the following 
correlations : 



(%„(*)%„,(*')) = 2DL6 n ^ m S(t - t'). 



(11) 



Again by decomposing Fourier modes rj qn = £ gri + iQ qn in 
their real and imaginary parts, one obtains, in addition 
to relations = £ q _ n , C?„ = ~~ Cq- n > conditions on the 
correlations 



(£*»(*)&» (*')} = DLS ntm S(t - t'), 
(«0) =DL5 n , m 5(t-t'), 



£„(*)<,„.(*')> =0 



(12) 

(13) 
(14) 



for n > 0, 
gets Cqo 



m > 0. 
and 



For the Fourier component n = 0, one 



L(t)L(t')) =2DL8(t-t') 



(15) 



Using Eq. ([)]) in Eq. (|l]), an infinite system of coupled 
Langevin equations is obtained : 



dh qn (t) 
dt 



The spectral approximation now amounts to project the 
above infinite system on the space of periodic functions 
of period L with a finite number of Fourier modes h Qrl 
{\ln\ < Qn/2)- All equations retain their original form 
with the proviso that infinite sums X)^L-oo are now re ~ 
placed by finite ones Y^„=-N/2- ^ ms P roce dure thus as- 
sumes that hg n = for any n > N/2, and the original 
continuum equation is then reduced to a set of N+ 1 real 
Langevin equations. The first one 



dh, 



A N/2 

dt cL+ iX^ 

n— 1 



[-1 



(17) 



governs the temporal evolution of the zeroth mode h qo 
and depends on h qn for < n < N/2. The other N real 
equations are independent of h qo and describe the time 
evolution of a qi , /3 qi ,..., a qN/2 ,0, 



1" 

<?N/2 



da„ 



dt 



d(3 q 



dt 



F„. 



G„ 



a, 



a, (3 



(18) 



(19) 



where 



F„ 



a, ft 



-V&q-n 



(20) 



N/2 



— y 

2L ^ 



Qm Qn 



771,771'=— N/2 



G 



qn 



a,/3 = -vq 2 n l3 qn 



2L 



N/2 

E 

m,m' = -N/2 



QmQw 



®q™Pq m > + Pqm^q m > 



(21) 



^n,m-\-m' • 



= cL6 nfi - vq\ h qn (t) 



A 

2L 



m.m' — — 



Qm Qm' hq m {t)h qm , (t)S n ,m+m' + Vq n (*) 



The philosophy underlying this regularization is akin 
a renormalization group scheme in momentum space, 
where high wavevectors are integrated out above a cut- 
off Qn/2 m momentum space. This approximation is an 
alternative - albeit not equivalent - method to discretize, 
in real space, the width L into N+l independent points 
separated by a distance a = L/N. For the sake of sim- 
plicity, the number N is hereafter assumed to be a power 
of two. 

For the continuum linear counterpart - the EW equa- 
tion - no approximations are involved in the framework of 
this spectral method, as equations for modes q n > q N / 2 
(16)are simply discarded. On the other hand, this approach 
is shown to be far more useful than the real space Eqs. (||) 
and (|j), in the non- linear KPZ case, since it preserves, 



unlike both approximations given in Eqs. (|) and 
some basic properties of the original KPZ continuum 
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equation as shown below. Indeed let us recall that, apart 
from a normalization factor, the steady state distribution 
of continuum KPZ is given by 



P s [h] ~ exp 



v 

2D 



L/2 



-L/2 



dx(d x h)' 



(22) 



The distribution of modes \q n \ < qpj/2 m momentum 
space thus reads : 



length scales a s > a can be simply obtained by cutting 
out modes with wavenumber larger than qN s /2 (N s — 
Na/a s ). From Eq. (p5|), it is straightforward to get 
the steady probability of these coarse-grained surfaces, 
that is a steady probability of the same form as equa- 
tion (HU) with N — N s and the same £ (L = N s a s ). 
The coupled Langevin equations and (|l9|), ensure 
that, under coarse-graining, the exact steady-state is re- 
covered. 



P,[h, 



exp 



N/2 
n=-N/2 



(23) 



III. EVOLUTION EQUATIONS FOR 
CORRELATION FUNCTIONS. 



The zeroth mode h qo does not contribute to ( p3[) s ince the 
average value of h(x,t) does not appear in (E2h, mean- 
ing that the surface always grows and its average never 
settles to a steady value unlike surface gradients. Wc 
note that both the original distribution Eq. ( p2[ ) and its 
spectral approximation Eq. (p3| ) are independent of A. As 
already remarked in Ref. |10||, this property is not satis- 
fied by the finite difference approximation Eq. (3). This 
inconvenient point is partly solved by the modification 
given in Eq. (||) which leads to the correct steady state 
distribution 



P.[hi 



IN 



exp 



N 
\ ^ 



h 3 ) 2 



(24) 



However it has been shown in Ref. [|12|, that even the 
approximation given in Eq. (^), fails in the presence of 
coarse-graining. This means that if fluctuating interfaces, 
obtained by the numerical generation of a real space dis- 
cretized KPZ equation at scale a, are smoothed out up to 
a scale a s > a, these coarse-grained surfaces cannot be 
described by a renormalizcd KPZ equation at a s . The 
origin of this problem can be traced back to the form 
of the steady state distribution (|24|) within a real space 
scheme. On the contrary, our spectral discretization has 
been devised in such a way to avoid this drawback (see 
below), and can then act as a safe starting point for the 
reconstruction procedure. 

On remembering that there are only N independent 
modes, the Fokker-Planck equation governing the evolu- 
tion of the probability distribution P[a,j3,t] associated 
with the Langevin equations (p~8|)-(p~9|) , reads Q: 



dP 
~dt 



N/2 

D 



d(F qn P) d{G qn P) 



da a 



d(3 Q 



DL 



d 2 P d 2 P . 



One may then check that the steady solution ( |23|) sat- 
isfies Eq. (^5|). It is known || that, when such 
a steady probability exists, all time dependent distri- 
butions asymptotically converge towards it. The dis- 
cretized KPZ thus preserves the particular symmetry 
of the continuum KPZ equation. A further advantage 
of this discretization, is that surfaces coarse-grained at 



Next we address the time dependent distribution 



P[a,f3,t] appearing in equation (|25j). We first exhibit 
an evolution equation for the ensemble average (3^) (t), 

and (d 2 gn )(t) 



(26) 



(a 2 q J (t) = / VaV/3 a 2 Qn P[a,0,t] 



Pl)(t)= I VaV/3 01 P[a,p,t], 



32 



where 



N/2 

VaVp = Y[ da qj d$ qj . 

3=1 



(27) 



(28) 



Assume that the probability density P[a,f3,t] goes ex- 
ponentially to zero for a, (3 going to 00. If the Fokker 
Planck equation ( p5| ) is multiplied by a 2 n and then inte- 
grated over all variables, one obtains, after an integration 
by parts : 



dt 



-I 



2/ VaD[3 a qn F qn [a,/3]P[a,P,t] + DL, (29) 



and similarly for p 2 (t), 



d{Pl){t) 



= 2/ VaV[3 p qn G qn [a,(3\P[a,f3,t] + DL. (30) 



dt 



From the expressions (|2l|)-(p2j), we then get for n > 



<»)i<|SKl" 



2v<ll(\h qn \ 2 )-\V qn +2DL, (31) 



where 



V„. = j ^2 q m q m 'TZ 

m, m' = -N/2 



N/2 



q " L 



h— q n h qrn h qrn , y ^n,m+m'? (32) 



and 1Z indicates that only the real part is considered. 
Equation ( |3l| ) implies that 
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-g 2 (t) = -2vg A {t) - XK(t) + 2DQ 2l 



where 



N/2 



92p( 



(t) = i J2 £ P (\KM\ 



(33) 



(34) 



n=-N/2 



for p = 1, 2, . . . and 



Jf(t) 



1 

ft 

N/2 



N/2 



^2 ^llrnQm' 
n,m,m f — — N/2 

h- qn (t)h qm {t)h q (t) 



(35) 



Q 2 



2 _ 4 2 7V(iV + l)(2iV+l) 



E2 * 



-—N/2 



I? 



(36) 



In Appendix [A] is shown that the term -fT(i) in Eq. j3^ ) 
actually vanishes identically. Hence one is left with 



dt 



g 2 (t) = -2vg 4 (t)+2DQ 2 . 



(37) 



Therefore only and D explicitly appear in the above 
equation. The absence of the parameter A in Eq. d37j ) 
is reminiscent of an analogous phenomenon appearing in 
the steady probability distribution (|2^) in 1 + 1 dimen- 
sional case. It is worth stressing that, although the A 
term does not explicitly appear in Eq. (|37j), it is never- 
theless implicitly present through the evolution of 54 (t). 

Another equation can be obtained by averaging 
Eq. © 



where 



j t 9o{t) = cL+^g 2 (t). 



ga(t) = (h go (t)) , 



(38) 



(39) 



This latter relation does not explicitly involve either the 
noise or the diffusion term. 



IV. THE PSEUDO-SPECTRAL METHOD FOR 
THE KPZ EQUATION 

A. Numerical procedures 

In order to compare rough surfaces numerically gener- 
ated by various spatial discretizations of the KPZ equa- 
tion, we introduce a one step Euler Scheme in time for 
the temporal discretization. This simple algorithm is 
used since, for stochastic equations, it is known that a 
naive application of a two- step s method can result in less 
computational efficiency fll9|] . In order to speed up the 



evolution of equations (|18|) and (|19|), a pseudo-spectral 
method is used at each time step of the Euler scheme. 
This method efficiently computes the quantity 



X_ 

2L 



N/2 

^ qmq m 'h qm (t)h qm ,(t)Sr h m+m', (40) 
m,m' = -N/2 



whose real and imaginary part are the non-linear con- 
tributions of the N + 1 Langevin equations ([l8|) and 
(|l9|). In this procedure, quantity Xq n can be evaluated 
without explicitly performing the double sum appearing 
in Eq. ( f4(i| ) for each n. First, the Fourier modes of the 
surface derivative d x h(x,t) are obtained by a simple al- 
gebraic multiplications iq n h qn for < n < N/2. One 
then returns to real space 



N/2 



d x h(x,t) = — ^ i( lnh qrl (t) e 



(41) 



-N/2 



to obtain the gradient at given spatial points. The 
computation of ^(d x h) 2 is then straightforward at these 
points. One then exploits again a Fourier transform to 
go back to spectral space using the values of the nonlin- 
ear terms at these collocation points. In order to prevent 
the aliasing problem |Q, we suitably choose the 

range of the collocation points and the number of Fourier 
modes to be used, in such a way that the procedure pro- 
vides the exact \ qn for < n < N/2 of Eq. (f40|). The 
modes external to the range —N/2 < n < N/2ca,ii be 
dropped since they do not enter in the actual computa- 
tions. This well-known technique (dealiasing |l4|]), albeit 
seemingly more complicated, turns out to be much more 
efficient of a brute force computation of Eq. (Eoh. 



B. Comparison with real space discretization 

We now compare the performance of the spectral dis- 
cretization with the one based on the standard Eq. (j^) 
and modified Eq. (^) real space discretization and with 
the corresponding analytical results obtained for the con- 
tinuum KPZ equation. The surface is grown via the 
pseudo-spectral KPZ equations with a regularization at 
scale a. In these simulations we have always used an Eu- 
ler time step in the range 10 -2 — 10~ 3 which is sufficiently 
small to avoid any numerical instability up to the sizes 
and for the A considered in here. 

In order to compare the performances of the pseudo- 
spectral method with respect to the real space discretiza- 
tion schemes, we apply a test akin to the one carried of 
in Ref. JO]]. We generate a steady state KPZ surface in 
(1 + 1) dimension, with lattice spacing a = 1, and pa- 
rameters v = 1, A = 3, D = 1, c = 0, by using (i) the 
standard real space discretization (|3 | ), (ii) the real space 
discretization ^ introduced in Ref. |11] , (iii) our pseudo- 
spectral discretization. The steady state roughness W(L) 



5 



obtained by the three methods for various sizes L is then 
compared with the exact value 



D 



w ^ = \lw„ Ll/2 > 



(42) 



of the continuum KPZ equation ^J. Fig. |T| shows that, 
unlike the standard real space representation (Q) which 
underestimates the ratio v/D (= 1 in the present case), 
both the modified real space representation (||) and the 
pseudospectral representation yield, on average, the cor- 
rect ratio. This is no surprise since we have previously 
shown that the pseudo-spectral representation correctly 
accounts for the steady-state distribution properties, a 
feature not shared by of the standard real space repre- 
sentation (|) @. 

Figure || depicts the behavior of gi (t) for a surface of 
a size L = 256, flat at t = 0, and average over TV = 500 
independent growths. The curve has a gradual increase 
(starting from zero) until it saturates after a character- 
istic time tf^ ~ 30000 x 10 -3 This time, which depends 
upon the size of the system, represents the cross-over af- 
ter which (i) all length scales have saturated, (ii) the ve- 
locity of the average height becomes constant and (iii) the 
roughness levels off. From renormalization group theory, 
one expects that i^ ~ L z where the dynamical exponent 
z is equal to 3/2 for 1 + 1 dimension 

In the spirit of renormalization group theory, let us 
separate the dynamics of modes of wavelength shorter 
than a s = ba with b = 2 l with modes of wavelength 
longer than a s . This can be done using the following 
decomposition of gi p (t) ( N s = N/b) 



92p(t) 



<$(*) + <(«), 



with 



sS?(«) 



N B /2 

E n 

n=-N B /2 



2p 



IV (t)| s 



(43) 



(44) 



representing the "slow" part, and 



R 



-N./2-1 



n=-JV/2 
N/2 

2 E #(fth.(*)l a )> 



n=N s /2+l 



indicating the "fast" part. In Fig. ^, the behavior of 



.92 (i) , an( i R^'it) are plotted for a scaling b = 2. 

Before a characteristic time ( t\ <~ 350 x 10~ 3 for 
the parameters chosen in the pseudo-spectral KPZ equa- 
tions), shortwave modes N s /2 + l < \n\ < N/2 contained 
in R 2 (t) are evolving much faster than longwave modes 
< \n\ < N s /2. Similar features occur for higher val- 
ues of b = 2 l thus defining a sequence of characteristic 



? (6) 



times t\ < t% < t% < 



1300 x 10~ 3 and t% 



For instance it is found that 
- 3000 x 10~ 3 . The above re- 



marks are clearly important when defining the dynamics 
of the coarse-grained surface obtained by eliminating the 
modes of wavelength shorter than scale a s = ba = 2 l a. 
This surface is characterized by the same average height 
= go for any b, and g^ (t) playing the role of 52 (t) ■ 
By rewriting Eq. (j^) as 



dt 



1+ 2L R 



(b) 



L + ±gi b \t), 



(46) 



it is clear that the coarse-grained surface with N s = N/b 
modes will satisfy a KPZ equation with a renormalized 
c s and identical A s = A only when B^' has reached 
saturation. On starting from a flat interface, this hap- 
pens whenever the Fourier modes N s /2 < n < N/2 and 



—N/2 < n < -N s /2 are thermalized, i.e. for t > tf. 
Similarly, upon summing Eq. ([ ' 
n = 1, ... , N s /2, one finds that, 



over the slow modes 



dg ( 2 b \t) 
dt 



where 



= -2 V g ( i\t)-\ 



N B /2 



'hi 



2DQ$\ 



(47) 



-N./2 



Qi b) = 



N s /2 

E & 

n=-N 3 /2 



(48) 



For t > t\ , it is assumed that dynamics of fast modes is 
slaved to the one of slow modes, in such a way that the 
term in A can be written as 



N 3 /2 

' E 

n=~N 3 /2 



<l,< 



= 2Avg?\t) 



2ADQ%\ 



(49) 



(45) |. 



thus satisfying a coarse-grained KPZ equation with 
renormalized parameters v s — v + Av and D s = D + AD. 
For the EW universality class, it is easy to show that, us- 
ing Eq. (46) and Eq. (f47j), all parameters do not renor- 
malize, as expected from renormalization group theory 



Our purpose here is to compute the KPZ renormalized 
parameters directly from the experimental observations 
of the surface growth. According to the previous dis- 
cussion, we should begin collecting the data after a time 
t > tf (starting from an initially flat surface) in order 
to best characterize the dynamics of a surface coarse- 
grained by a factor b — 2 l . Moreover the most prominent 
evolution of this surface occurs during the time interval 
when the length scales between 2a s and a s = ba are not 
thermalized i.e. during the interval [tf, if. J which is then 
the optimal period to characterize its dynamics. 
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V. THE RECONSTRUCTION PROCEDURE 
BASED ON THE SPECTRAL APPROXIMATION. 

Based upon the above spectral approximation, we now 
introduce a method to identify, at a given length scale, 
an optimal set of KPZ effective coupling parameters c, v, 
A and D, from "experimental" snapshots of interface pro- 
files. The "experimental" data may be generated by nu- 
merical simulations of the KPZ equation itself (see tests 
in section VI), by numerical m icros copic models emu- 
lating surface processes (Section VII ) or associated with 
real interface growth. For the sake of simplicity, mea- 
surements are always assumed free of observational noise. 
This approach thus constitutes a typical inverse problem 
for an infinite (or finite, but large, in the discretized ap- 
proximation) dimensional system with a finite number of 
parameters to identify. Although in this work we focus on 
the KPZ universality class, our method has a more gen- 
eral validity, and may be extended, with slight changes, 
to other universality classes. 

A first reconstruction of KPZ dynamics was attempted 
by Lam and Sanders ]lfi|] . These authors used equa- 
tion (|^) and worked in real space using experimental 
heights h° bs (t) measured at N points Xj — ja with 
j = 1,...N. For the reconstruction, they performed a 
least-squares calculation directly on the Langevin equa- 
tions to compute the parameters c, v, A. The noise term 
D was eventually obtained as a by product of the pre- 
vious calculation. In Ref. |0|, the difficulties associated 
with this approach have already been discussed. In the 
present work, the analysis is based on the philosophy ex- 
plained in previous Sections and equations (^7|), respec- 
tively (|38|), are used to identify through a least-squares 
procedure the coefficients v, D, respectively c, A. Be- 
sides the fundamental theoretical features discussed in 
Sec. [n], several reasons may be invoked in favor of our 
approach. First, this method is not directly based on the 
primitive stochastic equations but uses the deterministic 
equations introduced in Section III which govern the en- 



semble average of correlation functions. Standard identi- 
fication algorithms (e.g. the least-squares method) which 
are well suited and have been widely tested for deter- 
ministic equations, are then expected to be more reliable 
under such conditions. Second, as functions g2p(t) are 
already averaged quantities, a smaller number of realiza- 
tions is expected to be required due to the self-averaging 
property. Finally, dynamical noise is directly introduced 
in our reconstruction algorithm, unlike in the one of 
Ref. @. This seems a natural requirement since noise 
is an intrinsic parameter of the interface evolution. 

In order to get ensemble averages of spatial correla- 
tions, Af growths starting from the same surface, e.g. 
a flat surface, have been carried out. This provides Af 
distinct realizations of the same stochastic process. For 
a given realization, the experimental surface is observed 
at time (tk = kAt, k = 1,2, ..,pM) where At is the mea- 
surement sampling time. This procedure may be easily 



performed in a real experiment and leads to quantities 
9o(t)i 92(t) and (t) which are linked, at the scale con- 
sidered, to the average height, the average of the square 
of the first or second surface derivatives of the smoothed 
surface. 

Let us now integrate equation (^) during p sampling 
times At : 



g 2 {T k+ i) -92{T k ) 



-2v- 



1 



Tk+i — Tk J T 



gi (t)dt + 2DQ 2 , (50) 



where time Tk — ftp At, (k = 1, . . . , M). Similarly from 
equation (F58h one gets 



go (T t 



k+l) 



9o{Tk 



Tk+i — Tk 



1 



2 Tk+\ — Tk 



92(t)dt. (51) 



If At is smaller than the characteristic time of the dynam- 
ics, one may approximate the time integral in Eqs. j50| ) 
and @, as an average over the p intermediate sampling 
times, thereby obtaining M — 1 linear constraints on the 
parameters v, D and c, A. A simple least-squares calcu- 
lation then determines v, D from Eq. ( |5fj| ) and c, A from 
Eq. P). 



VI. RESULTS ON THE KPZ RECONSTRUCTION 

In order to test our reconstruction method, we use the 
above procedure on rough surfaces generated through a 
discretized KPZ equation with known coupling param- 
eters (y — 1,A = 3,D = 1, c = 0). Note that, in or- 
der to compute the ensemble average values go{t), g 2 {t) 
and g±(t), we use Af = 500 samples to obtain a con- 
vergent value. This is important since statistical fluc- 
tuations may be large enough to induce a measurement 
error exceeding the non-linear contributions, and under 
such circumstances the identification would fail. Further- 
more each reconstruction is repeated for a small number 
of independent configurations (typically 5) in order to get 
an estimate of the error bars associated with the recon- 
structed parameters. 

In the absence of coarse-graining, data are produced 
by the discretized equations (|l7|)-(|lS|)- ( |l9| ) where the 
minimum length scale a = L/N is introduced by the 
spectral approximation. Since equations used in the re- 
construction procedure are exactly identical to the ones 
generating the time series, this provides a first stringent 
test on the validity of our method. For the original data, 
the most efficient interval choice to perform the identifi- 
cation is expected to be [0, if]. Indeed, after t\ : all length 
scales between a and 2a are thermalized and the contri- 
bution to the variation of steady velocity, stemming from 
the non-linear term, becomes less and less important (see 
equation (^6j)). In this instance, the inversion technique 
may run into difficulties to compute c and A since statisti- 
cal fluctuations in the computation of ensemble averages 
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should not be greater than the value of this nonlinear 
contributions, as mentioned above. 

Such a procedure may be iterated for coarse-grained 
surfaces at a s — ba (b = 2 l = 2, 4, . . .) which are still 
assumed to be governed by KPZ-like equations. The re- 
construction of the renormalizcd equation should then 
be performed with data taken after tf since the renor- 
malized equation is not expected to be valid at earlier 
times. Once again, because of statistical fluctuations, we 
expect the most efficient interval choice for the identi- 
fication at length scale a s to be [if,if +1 ]. Clearly the 
reconstruction becomes more and more difficult to im- 
plement as b increases, since time intervals [tf, tf +1 J and 
required statistics will correspondingly increase. Fig.||, [|, 
H depict the results for the three parameters u s , A s and 
D s as a function of the scaling ratio b for various lattice 
sizes (error bars are of the same order of magnitude of 
the symbol sizes and therefore not shown). The opti- 
mal values is then expected to correspond to the L — > oo 
limit. For the scale a, one obtains the correct values 
which shows the capability of the method to identify the 
correct parameter. For the coarse-grain case the extrap- 
olated values of A, as well as the ratio v/D, appear to 
be independent of b (as they should) while parameters v 
and D are renormalized S. 

VII. RESULTS ON GROWTH MODELS 

The ultimate goal for our method is to be applied to 
real experimental surfaces. In this section, as an inter- 
mediate step, the reconstruction technique is tested on 
data produced through numerical microscopies models. 
Specifically we consider two typical models: (i) a random 
deposition model with surface diffusion (RDSD) which is 
described by the EW linear continuum theory (A = ); 
(ii) a particular solid-on-solid model called single-stepl 
(SSI) which is expected to belong to the KPZ universal- 
ity class. This latter model can be mapped onto an Ising 
model pl| ] thus providing an analytical value of A and 
hence a further stringent test for our method. 

A microscopic growth model is composed of three main 
ingredients. (LI) A probabilistic law, independent of the 
surface dynamics itself, which describes the flux of parti- 
cles directed towards the surface; (L2) a - deterministic 
or probabilistic - rule which determines whether a given 
particle directed towards the specific site s is effectively 
deposited (s is an active site) or simply discarded(s is not 
active); (L3) a prescription yielding the displacement of 
the particle or the rearrangement of the surface, after 
that the deposition has become effective. 

Law (LI), characterizing the particle flux, defines the 
mean time or the time scale necessary to a particle to fall 
towards - but not necessarily deposited on - the surface. 
In principle, this law may be given by a probability vary- 
ing with site locations and may also be of intermittent 
nature. In this work, however, we consider the simplest 



case of uniform flux both in space and time. During 
each time interval St, which defines the time scale of the 
process, a unique particle falls on the surface with prob- 
ability one and selects a given site s with equiprobability 
1/N. 

In order to compute the parameters, we use as space 
unit the unit cell of the microscopic model a = 1 {L = N 
where N is the number of sites). The unit of time is 
defined in such a way that St = 1/N. In the following, 
we use sizes up to L = 8192 with J\f = 50 independent 
growths. 

A. The RDSD model. 

In the RDSD, a particle directed towards the specific 
site s is always deposited (L2). This means that, one 
layer is deposited on the surface during a unit of time. 
Law (L3) can be phrased as follows. Upon reaching the 
surface, the particle falling towards a specific site s, com- 
pares the heights of the nearby neighbors and sticks to 
the one of lowest height unless the original site is a local 
minimum (in that case it does not move). 

We have generated a RDSD starting from a flat surface 
and performed coarse graining up to b = 16. The sam- 
pling measurement time At is taken to be larger then the 
time unit. We have consistently found A ~ and c ~ 1, 
as expected. Moreover we find a value for v and D which, 
as b increases, slowly converge to v ~ 0.8 and D ~ 0.5 
(we recall that c — 1 and D = 0.5 for a simple random 
deposition model with our time units). 

B. The SSI model. 

For the SSI model, an "active site" is defined j2^] as a 
site which is a local minimum for nearby sites, i.e. such 
that hi < hi-x and hi < hi+i (L2) . The (L3) rule for 
the SSI model imposes that two particles are deposited 
on the active site. Our algorithm has been devised to 
cover non-steady-state conditions. Since this situation is 
hardly discussed in the literature j|] , j2^] , in Appendix [B| 
the way of including the time dependence in this model in 
an efficient way, is reported. The sampling measurement 
time At is taken to be the unit of time. Note that, unlike 
the RDSD model, the time unit does not correspond to 
the effective deposition time of one layer. 

We have used a "teeth-like" initial surface with odd 
and even sites having height and 1 respectively (hence 
there are N/2 active site at this stage) |22j. Reported in 
Fig. |?] is the result for the ratio v s j D s for various scaling 
factors b. As b increases, v s /D s tends to a constant char- 
acteristic of the KPZ phenomenology for b > 4. This is 
confirmed by Fig. || which depicts the results for A s for 
b > 4, displaying a tendency for X s towards —2. 

A word of caution is in order here. As we discussed, 
the value of A is exactly known through a mapping onto 
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an Ising model p3]. The predicted value is A = — 1 when 
computed with a time unit yielding c = 1 at stationarity. 
This result is consistent with ours since it corresponds 
to a time unit which is half of the one we have exploited 
{c s ~0.5). 



VIII. CONCLUSION. 

In this work we have proposed a new method to extract 
the coupling parameters of the KPZ equation from exper- 
imental snapshots of successive interface profiles. This 
method hinges on two main ingredients. First a pseu- 
dospectral scheme is used to simulate the KPZ equation, 
and this scheme can be reckoned as an improved dis- 
cretization with respect to the standard real space finite 
difference ones. As a matter of fact, it preserves both the 
correct steady state distribution and the coarse-graining 
properties of the continuum corresponding equation. Sec- 
ond, our reconstruction algorithm is based on the time 
evolution of correlation functions. These functions do 
satisfy a deterministic evolution equation entitling the 
use of standard least-square procedure to identify the 
coupling parameters. This second ingredient parallels the 
analogous one introduced in Ref. [ fl2| which was however 
based on a real space representation. 

We have first tested the overall procedure on numer- 
ically generated KPZ profiles with known coupling pa- 
rameters. Our scheme is capable to reproduce not only 
the correct parameters in the absence of coarse-graining 
but unlike the previous attempt [ fl2| , it also provides con- 
sistent and robust results for coarse-grained surfaces. 

Next we have applied our algorithm to microscopic 
models which more closely mimic experimental situa- 
tions. In such smoothing procedure is unavoid- 
able in order to describe the surface in terms of a con- 
tinuum evolution equation, and it is then a vital require- 
ment to use an efficient and reliable method under such 
conditions. Again we were able to reproduce few known 
analytical results for these microscopic growth models. 
Furthermore, some additional estimates of other param- 
eters have also been given. 

We remark that our method is of general applicability. 
For instance, an extension to two-dimensional surfaces 
is not expected to present major difficulties. Similarly 
it could be applied to determine stochastic equations 
emulating coarse grained equations of the Kuramoto- 
Shivashinski type where the universality class is still an 
open question |23|. 
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APPENDIX A: PROOF THAT K(T) = 

We prove here that the quantity K(t) appearing in 
Eq. d33| ) is actually zero. The proof is patterned after a 
similar proof used to show that the steady state proba- 
bility d23| ) is independent of A for a (1 + 1) KPZ equation. 
Indeed by a reflection q„ — > —q n we obtain 



N/2 



m = j2 E 

k((v(«)V (t) 



(Al) 



-—N/2 



So 



n+m+m' ■ 



The above quantity is invariant under a cyclic permuta- 
tion of the indices n, to, to'. Therefore it can be rewritten 
as: 



1 



1 



N/2 

E ^qnq m q m '[qn + q m + q m >\U( (h qn (t)h qm (t)i 

n,m : m' — — N/2 



K(t) 

which obviously vanishes. 



APPENDIX B: TIME EVOLUTION OF THE SSI 
DEPOSITION MODEL: AN EFFICIENT 
ALGORITHM 

A naive simulation of the SSI model would proceed as 
follows. During each time interval 5t, a unique particle 
is dropped on the surface, and one checks whether it is 
falling on an active or inactive site according to (L2). If 
deposition is attempted on an active site, according to 
law (L3) for SSI, time is incremented and the particle is 
deposited. On the contrary, the particle directed towards 
an inactive site is discarded but the time is nonetheless in- 
cremented. This procedure is however very time consum- 
ing. Efficient algorithms, generating equilibrium surfaces 
in a fast way, simply consider active sites and do not take 
time into account. Such algorithms cannot be used here 
since (i) we have not reached equilibrium and (ii) we need 
to quantify the interface effective time evolution. There- 
fore we next implement an improved algorithm providing 
the time evolution as well. 

At time t = kSt, let us assume we know the number 
N a (t) of active sites of the surface and their respective 
positions. We then select, with equal probability, one of 
such active sites, we deposit a particle on it and then 
perform the re-ordering of the surface and the updating 
of the active sites. The important question is how long 
we have to wait to see the deposition event to occur. 
The probability of deposition between t and t + St is 
a = N a (t)/N and the probability that k — t/St time 
steps elapse before a deposition occurs is a[l — a] 1 . 
This law of probability of time intervals is numerically 
generated as follows. Intervals [#0,6*1], 62], ■■■[Ok, dk+i] 
are defined in [0, 1] where 9q = 0, 0\ = a and 9k+i — 
6k+a[l— a] k , i.e. 9k = 1— [1— a] k . Assume that a number 
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(3 is chosen with an equiprobability in the interval [0,1]. 
If it lies in the interval [9k,0k+i], then the waiting time 
is given by (k + l)6t. 
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FIG. 1. i)(L) = y/l2Zw(L) as a function of the 
size L where W(L) stands for the steady state rough- 
ness computed using the standard real space discretization 
Eq.(^), the modified real space discretization Eq.(^), and the 
Pseudo-spectral discretization given in Eq. © and Eq.© 
for a one- dimensional KPZ equation. The dotted line cor- 
responds to the exact continuum value given by Eq. ( [42] ) 
(yjD = 1). Units are arbitrary. 
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FIG. 2. Temporal behavior of g2(t) for a KPZ growth on 
one dimensional substrate of size L — 256. The temporal 
axis t is written in terms of numbers of Euler time steps, here 
taken to be 10" 3 . 
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FIG. 3. Temporal behavior of functions gi (t), (t) and 
R ( 2 ] (t) (see text) for a KPZ growth on one dimensional sub- 
strate of size L = 256. The temporal axis t is written in terms 
of numbers of Euler time steps 10~ 3 . 



FIG. 5. The coupling parameter A s for various b and in- 
creasing lattice sizes 2048, 4096 and 8192 in the case of a nu- 
merically generated KPZ equation. The input value is A = 3. 
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FIG. 4. The coupling parameter v a for various 
coarse-graining levels b = 2 l and for increasing lattice sizes 
2048, 4096 and 8192 in the case of a numerically generated 
KPZ equation. The input value is v = 1. Error bars are of 
the order of the symbol sizes. 



FIG. 6. The coupling parameter D s for various b and for 
increasing lattice sizes 2048, 4096 and 8192 in the case of 
a numerically generated KPZ equation. The input value is 
D = 1. 
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FIG. 7. The ratio R = u s /D a for the SSI growth model. 



FIG. 9. The coupling parameter c 3 for the SSI growth 
model. 
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FIG. 8. The coupling parameter A s for the SSI growth 
model. 
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